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Converting  GPS  Coordinates  (cf)Xh)  to  Navigation 
Coordinates  ( ENU ) 


EXECUTIVE  SUMMARY 

In  many  applications  relevant  to  the  Australian  Defence  Force  (ADF)  it  is  necessary 
to  convert  the  Global  Positioning  System  (GPS)  coordinates  of  latitude,  longitude  and 
height  to  a  local  navigation  frame  with  coordinates,  east,  north  and  up.  For  example  when 
testing  navigation  instruments,  such  as  the  inertial  navigation  system  (INS),  it  is  often 
helpful  to  compare  these  measurements  with  those  obtained  from  an  independent  GPS 
receiver.  An  INS  records  the  east,  north  and  up  displacement  from  its  point  of  origin.  To 
compare  INS  and  GPS  measurements  we  need  to  transform  GPS  coordinates  to  navigation 
coordinates.  Furthermore,  east,  north,  up  coordinates  are  essential  in  determining  the  line 
of  sight  for  terrain  data  given  as  latitude,  longitude  and  height,  such  as  digital  terrain 
elevation  data  (DTED).  For  large  amounts  of  data,  e.g,  trial  data,  this  process  may  be 
very  computationally  intensive. 

Means  for  converting  GPS  data  to  navigation  frame  coordinates  already  exist.  How¬ 
ever,  the  method  presented  in  this  report  is  roughly  three  times  faster  than  a  commonly 
employed  one.  The  coordinate  transformation  routine  outlined  in  this  report  is  accurate 
to  within  10m  over  a  range  of  60km. 

The  coordinate  transformation  method  outlined  here  will  be  used  in  the  Navwar  simu¬ 
lation  package  currently  under  development  in  DSTO  Edinburgh.  It  will  have  the  effect  of 
significantly  reducing  the  processing  time  of  line  of  sight  calculations  in  jamming  models. 
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1  Introduction 

The  aim  of  this  report  is  to  provide  a  useful  and  efficient  means  of  using  Global 
Positioning  System  (GPS)  data  to  determine  local  range  and  bearing  measurements.  This 
is  particularly  pertinent  owing  to  the  ubiquity  of  GPS  receivers. 

The  output  of  a  GPS  receiver  is  latitude,  longitude  and  height  in  the  World  Geodetic 
System  1984  (WGS84)  coordinate  frame.  However  in  many  applications  we  are  interested 
in  the  range  and  bearing  of  these  coordinates  with  respect  to  a  reference  point,  for  ex¬ 
ample  the  range  and  bearing  of  an  aeroplane,  with  respect  to  the  control  tower.  When 
testing  other  navigation  instruments,  such  as  the  inertial  navigation  system  (INS),  it  is 
often  helpful  to  compare  these  measurements  with  those  obtained  from  an  independent 
GPS  receiver.  An  INS  records  the  east,  north  and  up  displacement  from  its  point  of  ori¬ 
gin.  To  compare  INS  and  GPS  measurements  we  need  to  transform  WGS84  coordinates 
to  navigation  coordinates.  Furthermore,  (east,  north,  up)  coordinates  are  essential  in  de¬ 
termining  the  fine  of  sight  for  terrain  data  given  as  latitude,  longitude  and  height,  such  as 
digital  terrain  elevation  data  (DTED). 

At  this  point  the  reader  may  ask,  “hasn’t  this  already  been  done?”  Yes,  it  has  been 
done  before,  for  example  GPSoft™  have  Matlab™  toolboxes  which  convert  GPS  coordi¬ 
nates  into  local  navigation  coordinates.  In  this  report  we  derive  a  method  which  does  this 
conversion  three  times  faster  than  GPSoft™.  This  difference  is  significant  when  dealing 
with  large  amounts  of  data. 


2  Converting  from  WGS84  to  navigation 

coordinates 


Navigation  coordinates  are  determined  by  the  fitting  of  a  tangent  plane  to  a  fixed  point 
on  the  surface  of  the  Earth.  The  e  axis  points  east,  the  n  axis  north,  and  u  axis  points 
perpendicular  to  the  tangent  plane  and  away  from  the  centre  of  the  Earth. 


Figure  1:  Navigation  and  ECEF  frames. 


What  follows  below  is  a  derivation  of  the  transformation  from  latitude,  longitude  and 
height  ( <j)\h )  to  east,  north  and  up  (ENU).  It  has  been  included  for  completeness  but  the 
reader  may  safely  jump  to  equation  (4)  if  they  are  only  interested  in  the  final  result. 
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The  transformation  from  4>Xh  to  ENU  is  a  three  stage  process: 

1.  Determine  latitude,  longitude  and  height  of  reference  point,  e.g.  location  of  radar. 
In  this  report  it  is  the  arbitrary  point  (<p,X,h). 

2.  Express  small  changes  in  latitude,  longitude  and  height  in  Earth  Centred  Earth 
Fixed  (ECEF)  coordinates. 

GPS  coordinates  can  be  converted  into  ECEF  coordinates  using  the  following  for- 


mulae: 

x  —  ( ^  +  h)  cos  <f)  cos  A 
y  =  (§  +  M  cos  ^6  sin  A 

2  =  ^  +  h)  sin  0 

(1) 

Where 

i  0 

X=  \J  1  -  e2  sin2  <p, 

a  and  e2  are  the  semi-major  axis  and  the  first  numerical  eccentricity  of  the  Earth 
respectively. 


To  convert  small  changes  in  latitude,  longitude  and  height  into  ECEF  coordinates  we 
need  to  Taylor  expand  equation  (1)  about  <f>  ->  <p  +  d<f>,  A  ->  A  +  dX,  and  h  h  +  dh. 

dx  =  (-acosAg0(i— e2)  -  hcos  Asimft)  d<f>  —  ( Q sin  ^cos ^  +  h sin  A  cos  d\ 

+  cos  cj>  cos  Xdh  +  ^acoscj)cos  A  (—2  —  7e2  +  9e2cos2<^>)  —  cos  A  cos  d$ >2 

+  /asm  Asin^>(i  e  )  ^  hsiuXsiruf)')  d<fidX  —  cos  A  sin  4>dhd4> 

+  ?_acos\cos<i>  _  i  ^  cog  ^  cog  ^  _  sjn  \  cos  (j)dhd\  +  0(dd3)  +  0(dhdd2) 

dy  =  jzgy-^yi-e1)  -^sinAsin^,)  tUf>  +  +  hcos  Xcoscf>)  dX 

4-  sin  <p  cos  Xdh  +  f|acos^sinA  (—2  —  7e2  +  9e2cos2  <fi)  —  sin  A  cos  (f^J  d<j>2 
+  (  — — °s  *  ^3  ^ 1  e  ^  —  h  cos  A  sin  0^  d(f>dX  —  sin  A  sin  4>dhdcf> 

+  (^■-  si2yCOS<*  -  \ h  sin  A  cos  d>)  dX2  +  cos  A  cos  cpdhdX  +  O(d0 3)  +  O(dhd02) 

dz  =  f  ^ . ^  jcos^  +  h  cos  4>)  )  d(f)  +  sin  <pdh  +  cos  (pdhdfb 

+  ^asin0  (—2  —  e2  +  9e2  cos2  <f>)  —  ^hsmcfij  d4> 2  +  0{dhd<p2), 

(2) 

where  d9  is  either  d(p  or  dX. 

3.  By  means  of  a  rotation,  displacements  in  ECEF  coordinates  are  transformed  to  ENU 
coordinates. 

The  ECEF  coordinates  (dx,  dy,  dz)  are  orientated  in  such  a  way  that  from  the  centre 
of  the  Earth  the  z  points  in  the  direction  of  true  north,  x  points  in  the  direction 
of  the  prime  meridian  and  the  direction  of  y  is  90°  from  the  prime  meridian,  see 
figure  1.  The  orientation  of  ENU  coordinates  is  determined  by  rotating  the  ECEF 
coordinates;  firstly  about  the  z  axis  by  A  degrees  and  then  the  new  y  axis  by  <fi 
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degrees, 


(  de  \  / 

'  —  sin  A 

cos  A 

0 

dn  = 

—  sin  cf)  cos  A 

—  sin  (f)  sin  A 

cos  <j> 

^  du  )  \ 

^  cos  4>  cos  A 

cos  4 >  sin  A 

sin</> 

(  dx  \ 

dy 

\dz  ) 


(3) 


Substituting  equation  (2)  into  equation  (3),  and  ignoring  terms  of  O(d03)  and 
0(dhd92)  and  higher,  we  get 


de  = 

(x  cos<t)dX  —  f— +  h^j  sin  4>d4>d\  +  cos  cf>dXdh 

dn  = 

(a(1x^  +  h^j  d(f)  +  |  a  cos  (f>  sin  (f>e2d(f)2  +  dh  d4> 

+ 

^  sin(/>cos  4>  ^  +  h  j  dX2 

(4) 

du  = 

dh  —  (l  —  f  e2  cos  4>  +  dcj)2 

- 

1  ^  a  cos2  <j>  hcos2^dX2 

Note: 

The  coordinates  (<f>,\,h)  are  ellipsoidal  (WGS84),  not  spheroidal  (geocentric).  The 
difference  is  most  easily  explained  by  figure  2. 


Figure  2:  Ellipsoidal  and  spheroidal  coordinates 


3  Matlab  Program:  dllh2denu 

Armed  with  equation  (4)  we  are  now  in  a  position  to  write  a  Matlab  program.  The 
code  for  our  program  dllh2denu.m  is  given  in  appendix  A.  The  routine  dllh2denu  is  simple 
to  use,  it  requires  two  inputs;  the  first  is  a  vector  indicating  the  location  of  the  reference 
point  (e.g.  radar)  in  latitude,  longitude  in  degrees  and  height  in  metres.  The  first  input 
is  a  1  x  3  vector  of  the  form 

IlhO  =  [0,  A,  h) . 
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The  second  input  is  a  n  x  3  matrix  indicating  the  location  of  the  points  of  interest  (e.g. 
the  location  of  the  aeroplane).  The  second  input  is  of  the  form 


4>  i  Ai  h\ 


llh  = 


4>n  An 


h 


71 


The  resulting  output  is  a  n  x  3  matrix  with  the  first,  second,  and  third  columns,  repre¬ 
senting  the  east,  north  and  up  displacement  respectively  in  metres 


denu  = 


e\  7i\  u\ 


€.n  Tln  Un 

As  an  example  suppose  the  location  of  the  reference  point  is 

[0  =  39°,  A  =  — 132°,  /i  =  0] 
and  three  points  of  interested  are  located  at 


A 

h 

39°  +  0.5° 

-132° 

0 

39°  +  0.5° 

-132°  +  0.5° 

0 

39°  +  0.5° 

-132° +  0.5° 

60000 

To  determine  the  location  of  these  points  in  ENU  coordinates  we  would  run  dllh2denu  in 
the  following  way: 


»  format  short  g 
»  IlhO  =  [39,  -132,  0] 

IlhO  = 

39  -132  0 

»  llh  =  [39  +  0.5,  -132,  0; 
39  +  0.5,  -132+0.5,  0; 

39  +  0.5,  -132+0.5,  60000] 

llh  = 


39.5 

-132 

0 

39.5 

-131.5 

0 

39.5 

-131.5 

60000 
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>>  denu  =  dllh2denu(llh0,llh) 
denu  = 

0  55510  -242.2 
43008  55629  -389.07 
43415  56153  59611 


4  Accuracy  of  Method 

The  displacements  as  determined  by  Eq.  (4)  are  approximations,  since  the  Taylor 
expansion  of  Eq.  (1)  was  truncated  at  second  order.  To  estimate  the  error  of  this  approx¬ 
imation  it  suffices  to  take  an  upper  estimate  on  the  effect  of  third  order  terms  O(d0 3)  and 
0(dhd62),  where  dd  is  either  d<j>  or  dX.  The  greatest  possible  third  order  effect  is  at  the 
equator.  Third  order  terms  at  this  point  are  of  the  order  add 3  and  dhdO2 .  If  we  limit  dO 
to  0.5°and  dh  to  60km  then  this  effect  is  less  than  7m.  Near  the  equator  an  angular  offset 
of  latitude  or  longitude  of  0.05°  equates  to  roughly  60km.  Hence  as  a  rule  of  thumb 

Equation  (4),  and  the  hence  the  program  in  appendix  A,  may  be  used  if  an 
error  of  less  than  10m  is  acceptable  and  the  range  of  measurements  is  less  then 
60km. 

Should  the  range  exceed  60km,  or  an  accuracy  much  less  than  10m  be  required,  it  is 
necessary  to  calculate  dx,dy,dz  using  the  method  employed  by  GPSoft™.  This  method 
firstly  converts  each  point  to  aq ,yi,Zi  using  Eq.  (1).  The  displacements  in  ENU 

coordinates  are  then  calculated  by  using  Eq.  (3),  and  the  fact  that  dxi  —  aq  —  xo,  where 
xq  is  the  location  of  the  reference  point  in  ECEF  coordinates. 

To  validate  our  code  let  us  compare  ENU  displacements  as  determined  by  GPSoft™’s 
llh2xyz.m  and  x yz2enu.m  (ENU qps0 ft)  with  those  determined  by  dllh2denu.m  (ENUsam). 
Taking  the  same  reference  and  data  points  as  before  we  find  that 


0,  A,  h 

ENUcpsoft 

ENUsam 

RMS  difference 

39°,  -132°,  0 

[0,  0,0] 

[0,  0,  0] 

0m 

39+0.5°,  -132°,  0 

[0, 

55509.42,-242.21] 

[0, 

55510.13,-242.20] 

0.7m 

39+0.5°,  -132  + 
0.5°,  0 

[43006.16, 

55627.52,-388.04] 

[43008.36, 

55629.06,-389.07] 

2.88m 

39+0.5°,  -132  + 
0.5°,  60000 

[43410.18, 

56152.22,59608.30] 

[43415.27, 

56152.66,59610.93] 

5.75m 

5  Efficiency  of  Method 


As  already  mentioned  in  the  introduction,  the  main  motivation  behind  this  work  is 
to  speed  up  the  processing  of  data.  In  this  section  a  simple  time  comparison  experiment 
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shows  that  data  can  he  processed  roughly  three  times  faster  using  the  method  outlined  in 
this  report. 

To  perform  this  comparison  fairly  the  programs  llhSxyz.m  and  xyzSenu.m  had  to  be 
modified  so  that  they  could  process  more  than  one  coordinate  at  a  time.  The  program 
comparison. m  was  used  to  compare  the  CPU  time  taken  to  calculate  ENU  using  these  two 
methods,  the  code  is  in  appendix  A.  The  results  are  shown  below. 

>>  comparison 

t_sat  =  0.935  +-  0.01354 

t_sam  =  0.344  +-  0.005164 

We  can  see  from  these  results  that  the  method  presented  in  this  report  is  roughly  three 
times  faster  than  that  used  by  GPSoft™1. 


6  Conclusion 

The  aim  of  this  note  is  to  report  work  in  DSTO  to  establish  an  efficient  method  by 
which  to  transform  GPS  coordinates  to  local  navigation  coordinates  east,  north  and  up. 
To  this  end  equation  (4)  was  derived  and  a  Matlab  program  dllhSdenu.m  was  created  to 
evaluate  it  for  any  number  of  data  points. 

Equation  (4),  and  consequently  the  program  dllh2denu.m ,  are  valid  only  if 

a)  The  range  is  less  than  60km,  and 

b)  an  error  of  a  few  meters  is  acceptable. 

The  accuracy  and  efficiency  was  compared  with  the  method  of  GPSoft™ and  a  simple  test 
showed  that  they  agreed  to  within  a  metre,  but  our  method  was  three  times  faster  than 
theirs. 
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Appendix  A:  Matlab  Code 


A.l  dllh2denu.m 

function  denu  =  dllh2denu(llh0 , llh) 

•/.•/.•/••/.•/.•/.•/.•/.’/.•/.CONSTANTS 
a  =  6378137; 
b  =  6356752.3142; 
e2  =  1  -  (b/a) ~2; 

•/.•/.’/.'/.■/.'/.•/.•/.•/.'/.Location  of  reference  point  in  radians 
phi  =  IlhO ( 1 ) *pi/ 180 ; 
lam  =  llh0(2)*pi/180; 
h  =  llh0(3) 

•/.'/.'/.•/.•/.'/.•/.•/.•/.'/.Location  of  data  points  in  radians 
dphi=  llh( : , l)*pi/180  -  phi; 
dlam=  llh( : ,2)*pi/180  -  lam; 
dh  =  llh( : ,3)  -  h; 

•/.•/.’/.•/.‘/.•/."/.•/.•/.•/.Some  useful  definitions 
tmpl  =  sqrt(l-e2*sin(phi)~2) ; 
cl  =  cos (lam) ; 
si  =  sin(lam); 
cp  =  cos (phi); 
sp  =  sin (phi); 

'/.•/.'/.'/.•/.•/.‘/.•/.•/.•/.Transformations 

de  =  (a/tmpl+h) *cp*dlam  -  (a*(l-e2)/(tmpl~3)+h)*sp.*dphi.*dlam  +cp. *dlam. *dh; 

dn  =  (a*  (l-e2)/tmpl,'3  +  h)*dphi  +  1 . 5*cp*sp*a*e2*dphi .  ~2  +  sp~2 .  *dh .  *dphi  +  ... 
0 . 5*sp*cp* (a/tmpl  +h)*dlam.~2; 

du  =  dh  -  0 . 5* (a-1 . 5*a*e2*cp"2+0 . 5*a*e2+h) *dphi . "2  -  ... 

0 .5*cp“ 2* (a/tmpl  -h)*dlam. ~2; 

denu  =  [de ,  dn ,  du] ; 
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A.  2  comparison,  m 


clear  all 

N=100000;  '/.  Number  data  points 
loops  =  10; 

mmmnrara*/.initiaiization 
phiO  =  39; 
lamO  =  -132; 
hO  =  0; 

orgllh  =  [phiO , lamO , 0] ; 

phi  =  zeros (N, 1) ; 

lam  =  zeros(N, 1) ; 

h  =  zeros(N, 1) ; 

enu  =  zeros (N, 3); 

xyz  =  zeros (N, 3) ; 

denu_sat  =  zeros(N,3); 

denu_sam  =  zeros(N,3); 

tl  =  zeros(loops,l) ; 

t2  =  zeros (loops, 1) ; 

mmm,/.mm,/.Coordinate  data 

for  i  =  1 : N ; 

phi(i)  =  phiO  +  0.5*i/N; 
lam(i)  =  lamO  +  0.5*i/N; 
h(i)  =  hO  +i; 
end 

llh  =  [phi,  lam,  h] ; 

Comparison  loop 

for  j  =  1: loops 
t  =  cputime; 

orgxyz  =  llh2xyzmodi (orgllh) ; 
xyz  =  llh2xyzmodi(llh) ; 
denu_sat  =  xyz2enumodi (xyz, orgxyz) ; 
tl(j)  =  cputime  -t; 
t  =  cputime; 

denu_sam  =  dllh2denu(orgllh, llh) ; 
t2(j)  =  cputime  -t; 
end 

xmxmxmxm  statistics 

tlm  =  mean(tl) ; 
tlstd  =  std(tl) ; 
t2m  =  mean(t2); 
t2std  =  std(t2); 

disp([’t_sat  =  » ,num2str(tlm) ,  »  +-  » ,  num2str(tlstd)] ) ; 
disp([’t_sam  =  * ,num2str(t2m) ,  ’  +-  » ,  num2str(t2std)] ) ; 
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Appendix  B:  Glossary,  Physical  Constants  and 
Mathematical  notation 


Symbol 

Value 

Meaning 

a 

6378137 

length  of  Earth’s  semi-major  axis  in  metres 

b 

6356752.3142 

length  of  Earth’s  semi-minor  axis  in  metres 

e2 

6.69437999013  x  10~3 

first  numerical  eccentricity 

4> 

latitude 

A 

longitude 

h 

height 

e 

displacement  in  the  east  direction  of  the  navigation  frame 

n 

displacement  in  the  north  direction  of  the  navigation  frame 

u 

up  displacement  in  the  navigation  frame 

X 

displacement  in  the  x  direction  of  the  ECEF  frame 

y 

displacement  in  the  y  direction  of  the  ECEF  frame 

z 

displacement  in  the  z  direction  of  the  ECEF  frame 
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